﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using SlimDX;
using System.Threading.Tasks;

namespace Blake.Fluid
{
    public class FluidSystem : PointSet
    {
        #region Fields

        float smoothing2;
        float poly6Kern;
        float spikyKern;
        float laplacianKernel;

        const float epsilon = 0.00001f;

        #endregion

        #region Properties

        public bool IsFountainOn { get; set; }

        #endregion

        #region Events

        #endregion

        #region Constructors

        public FluidSystem()
        {

        }

        #endregion

        #region Public Methods

        public override void Initialize(int numParticles)
        {
            base.Initialize(numParticles);

            particles.Clear();

            FluidSystemParams parameters = GetDefaultParameters();
            parameters.NumParticles = numParticles;

            Initialize(parameters);

        }

        public void Initialize(FluidSystemParams parameters)
        {
            //Clear buffers
            //Add buffer

            this.Parameters = parameters;

            ComputeKernels();

            Reset();

        }

        public override void Run()
        {
            // simulation scale (not Schutzstaffel)
            float ss = (float)(Parameters.ParticleDistance / Parameters.SimScale);

            if (Parameters.EmitterRate.X > 0 && (++frame) % (int)Parameters.EmitterRate.X == 0)
            {
                frame = 0;
                Emit(ss);
            }

            Grid_InsertParticles();

            var options = new ParallelOptions()
            {
                MaxDegreeOfParallelism = 8
            };

            //UpdateFixedParticles();

            //float minPressure = float.MaxValue;
            //float maxPressure = float.MinValue;

            Parallel.For(0, particles.Count, options, i =>
                //for (int i = 0; i < particles.Count; i++)
                {
                    SPH_ComputePressureGrid(i);
                    //float pressure = particles[i].Pressure;
                    //if (pressure < minPressure)
                    //    minPressure = pressure;
                    //if (pressure > maxPressure)
                    //    maxPressure = pressure;
                }
            );

            Parallel.For(0, particles.Count, options, i =>
                //for (int i = 0; i < particles.Count; i++)
                {
                    SPH_ComputeForceGridNC(i);
                }
            );
            Parallel.For(0, particles.Count, i =>
                {
                    Advance(i);
                }
            );
            time += DeltaTime;
        }

        public void CreateExample()
        {
            Reset();
            float sides = 40f;
            Vector3 volumeMin = new Vector3(-sides, -sides, -sides);
            Vector3 volumeMax = new Vector3(sides, sides, sides);

            //CreateWall(new Vector3(-4, -30, -8), new Vector3(-1.5f, 30, 8));
            //CreateWall(new Vector3(1.5f, -30, -8), new Vector3(4, 30, 8));

            CreateWall(new Vector3(-7, -30, -sides+8), new Vector3(-4f, 30, -sides+30));
            CreateWall(new Vector3(4f, -30, -sides+8), new Vector3(7, 30, -sides+30));

            ComputeKernels();
            Parameters.VolumeMin = volumeMin;
            Parameters.VolumeMax = volumeMax;
            //Parameters.InitMin = initMin;
            //Parameters.InitMax = initMax;
            //Parameters.PlaneGravityEnabled = false;
            Parameters.SimSize = Parameters.SimScale * (volumeMax.Z - volumeMin.Z);
            Parameters.ParticleDistance = (float)Math.Pow(Parameters.ParticleMass / Parameters.RestDensity, 1 / 3.0);

            float factor = 0.65f; //0.87f 
            float spacing = Parameters.ParticleDistance * factor / Parameters.SimScale;
            
            //var point = GetPoint(AddPoint());
            //var point2 = GetPoint(AddPoint());

            //point.Position = new Vector3(0, 0, 0);
            //point2.Position = new Vector3(spacing*0.28f, 0, 0);

            //Vector3 initMin = new Vector3(-10, 0, -10);
            //Vector3 initMax = new Vector3(-4, 0, 5);
            Vector3 initMin = new Vector3(-sides, 0, -sides);
            Vector3 initMax = new Vector3(sides, 0, -sides+20f);
            
            AddVolume(initMin, initMax, spacing);
            
            //initMin = new Vector3(5, 0, -sides);
            //initMax = new Vector3(sides, 0, 10);

            //AddVolume(initMin, initMax, spacing);

            float cellSize = Parameters.SmoothingRadius * 2.0f;

            Grid_Setup(volumeMin, volumeMax, Parameters.SimScale, cellSize, 1.0f);

            Grid_InsertParticles(); 
            UpdateFixedParticles();
        }

        public override void AddVolume(Vector3 min, Vector3 max, float spacing)
        {
            FluidParticle point;
            float dx = max.X - min.X;
            float dy = max.Y - min.Y;
            float dz = max.Z - min.Z;

            float randomDist = spacing * .1f;

            for (float z = min.Z; z <= max.Z; z += spacing)
            {
                for (float y = min.Y; y <= max.Y; y += spacing)
                {
                    for (float x = min.X; x <= max.X; x += spacing)
                    {
                        point = GetPoint(AddPointReuse());
                        point.Position = new Vector3(x + (float)random.NextDouble() * randomDist - randomDist / 2,
                                                     y,
                                                     z + (float)random.NextDouble() * randomDist - randomDist / 2);
                        point.Color = new Color3((x - min.X) / dx,
                                                 (y - min.Y) / dy,
                                                 (z - min.Z) / dz);
                        float zNorm = (z - min.Z) * (max.Z - min.Z);
                        point.Pressure = Parameters.RestDensity * (zNorm * 2);
                        //point.Pressure = (Parameters.RestDensity * (zNorm*2f + 1.0f)) * Parameters.InternalStiffness;
                    }
                }
            }
        }

        private void CreateWall(Vector3 min, Vector3 max)
        {
            Vector3 position = (max + min) / 2;
            RectangleParticle body = new RectangleParticle()
            {
                Position = position,
                Min = min - position,
                Max = max - position,
                IsFixed = true,
                Color = new Color3(1, 1, 1)
            };
            this.bodies.Add(body);
        }

        private void UpdateFixedParticles()
        {
            float factor = 1.0f; //0.87f 
            float spacing = Parameters.ParticleDistance * factor / Parameters.SimScale;
            spacing = factor * Parameters.SmoothingRadius / Parameters.SimScale;

            particles.RemoveAll(p => p.IsFixed == true);

            foreach (var body in bodies)
            {
                for (float x = body.MinPosition.X; x < body.MaxPosition.X; x += spacing)
                {
                    for (float z = body.MinPosition.Z; z < body.MaxPosition.Z; z += spacing)
                    {
                        AddFixedParticle(x, z);
                    }
                    //AddFixedParticle(x, body.MinPosition.Z);
                    //AddFixedParticle(x, body.MaxPosition.Z);
                }
                //for (float z = body.MinPosition.Z; z < body.MaxPosition.Z; z += spacing)
                //{
                //    AddFixedParticle(body.MinPosition.X + spacing, z);
                //    AddFixedParticle(body.MaxPosition.X - spacing, z);
                //}

            }
        }

        private void AddFixedParticle(float x, float z)
        {
            var point = GetPoint(AddPoint());
            point.Position = new Vector3(x, 0, z);
            point.Color = new Color3(0, 1, 0);
            point.IsFixed = true;
        }

        private void Advance(int particleIndex)
        {
            Vector3 min = Parameters.VolumeMin;
            Vector3 max = Parameters.VolumeMax;
            float SL = Parameters.Limit;
            float SL2 = SL * SL;

            float stiff = Parameters.ExternalStiffness;
            float damp = Parameters.ExternalDamping;
            float radius = Parameters.ParticleRadius;
            float ss = Parameters.SimScale;

            //foreach (var particle in particles)
            //Parallel.ForEach(particles, particle =>
            {
                var particle = particles[particleIndex];

                if (particle.IsFixed)
                {
                    return;
                }

                float adjustment, diff;
                Vector3 norm = new Vector3();
                particle.Acceleration = new Vector3();
                var acceleration = particle.SPHForce;
                acceleration *= Parameters.ParticleMass;

                float speed = acceleration.X * acceleration.X +
                              acceleration.Y * acceleration.Y +
                              acceleration.Z * acceleration.Z;
                if (speed > SL2)
                {
                    acceleration *= (float)(SL / Math.Sqrt(speed));
                }

                //fountain
                if (IsFountainOn &&
                    particle.Position.X > -4 && particle.Position.X < 4 &&
                    particle.Position.Z > -90 && particle.Position.Z < -10)
                {

                    var force = new Vector3(0, 0, 15.0f);
                    acceleration += force;
                }

                //Boundary conditions

                //Z-axis walls
                diff = 2 * radius - (particle.Position.Z - min.Z - (particle.Position.X - Parameters.VolumeMin.X) * Parameters.ZMin_Slope) * ss;
                if (diff > epsilon)
                {
                    norm = new Vector3(-Parameters.ZMin_Slope,
                                       0,
                                       1.0f - Parameters.ZMin_Slope);
                    adjustment = stiff * diff - damp * DotProduct(norm, particle.VelocityEvaluative);
                    Vector3 acc = norm * adjustment;
                    acceleration += acc;
                }

                diff = 2 * radius - (max.Z - particle.Position.Z) * ss;
                if (diff > epsilon)
                {
                    norm = new Vector3(0, 0, -1);
                    adjustment = stiff * diff - damp * DotProduct(norm, particle.VelocityEvaluative);
                    Vector3 acc = norm * adjustment;
                    acceleration += acc;
                }

                // X-axis walls
                if (!Parameters.WrapX)
                {
                    diff = 2 * radius - (particle.Position.X - min.X +
                           ((float)Math.Sin(time * 10.0f) - 1 + (particle.Position.Y * 0.025f) * 0.25f) *
                           Parameters.XminSin) * ss;

                    if (diff > epsilon)
                    {
                        norm = new Vector3(1, 0, 0);
                        adjustment = (Parameters.XminSin + 1) * stiff * diff - damp * DotProduct(norm, particle.VelocityEvaluative);
                        acceleration += norm * adjustment;
                    }

                    diff = 2 * radius - (max.X - particle.Position.X +
                           ((float)Math.Sin(time * 10.0f) - 1) * Parameters.XmaxSin) * ss;

                    if (diff > epsilon)
                    {
                        norm = new Vector3(-1, 0, 0);
                        adjustment = (Parameters.XmaxSin + 1) * stiff * diff - damp * DotProduct(norm, particle.VelocityEvaluative);
                        acceleration += norm * adjustment;
                    }
                }

                //Y-axis walls
                diff = 2 * radius - (particle.Position.Y - min.Y) * ss;
                if (diff > epsilon)
                {
                    norm = new Vector3(0, 1, 0);
                    adjustment = stiff * diff - damp * DotProduct(norm, particle.VelocityEvaluative);
                    acceleration += norm * adjustment;
                }

                diff = 2 * radius - (max.Y - particle.Position.Y) * ss;
                if (diff > epsilon)
                {
                    norm = new Vector3(0, -1, 0);
                    adjustment = stiff * diff - damp * DotProduct(norm, particle.VelocityEvaluative);
                    acceleration += norm * adjustment;
                }

                foreach (var body in bodies)
                {
                    CollideBodyParticle(body, particle);
                }

                //Gravity
                if (Parameters.PlaneGravityEnabled)
                {
                    acceleration += Parameters.PlaneGravityDirection;
                }
                if (Parameters.PointGravityEnabled)
                {
                    norm = particle.Position - Parameters.PointGravityPosition;
                    norm.Normalize();
                    norm *= Parameters.PointGravityStrength;
                    acceleration -= norm;
                }

                particle.Acceleration += acceleration;

                acceleration = particle.Acceleration;

                particle.Acceleration = acceleration;

                particle.IntegratePosition(DeltaTime);

                var velocityNext = particle.Velocity * DeltaTime / ss;

                // p(t+1) = p(t) + v(t+1/2) dt
                particle.Position += velocityNext;

                if (Parameters.ColorMode == ParticleColorMode.Velocity)
                {
                    adjustment = (Math.Abs(velocityNext.X) + Math.Abs(velocityNext.Y) +
                                 Math.Abs(velocityNext.Z)) * 20.0f;
                    if (adjustment > 1.0)
                        adjustment = 1.0f;
                    //adjustment = adjustment * 0.8f + 0.2f;
                    particle.Color = new Color3(adjustment, 0, 1 - adjustment);
                }
                else if (Parameters.ColorMode == ParticleColorMode.Pressure)
                {
                    float v = 0.5f + (particle.Pressure / 1500.0f);
                    if (v < 0.1)
                        v = 0.1f;
                    if (v > 1.0)
                        v = 1.0f;
                    particle.Color = new Color3(v, 0, 1 - v);
                }

                if (Parameters.WrapX)
                {
                    diff = particle.Position.X - (Parameters.VolumeMin.X + 2);
                    if (diff <= 0)
                    {
                        particle.Position.X = (Parameters.VolumeMax.X - 2) + diff * 2;
                        particle.Position.Z = 10;
                    }
                }

            } //foreach particle
            //);
        }

        #endregion

        #region Private Methods

        void CollideBodyParticle(RectangleParticle body, FluidParticle particle)
        {
            Vector3 norm;
            float adjustment;

            float stiff = Parameters.ExternalStiffness;
            float damp = Parameters.ExternalDamping;
            float ss = Parameters.SimScale;
            float radius = Parameters.ParticleRadius;
            float r2 = radius * 1.5f / ss;

            float diffXMin = 2 * radius - (body.MinPosition.X - particle.Position.X) * ss;
            float diffXMax = 2 * radius - (particle.Position.X - body.MaxPosition.X) * ss;
            float diffZMin = 2 * radius - (body.MinPosition.Z - particle.Position.Z) * ss;
            float diffZMax = 2 * radius - (particle.Position.Z - body.MaxPosition.Z) * ss;

            if (diffXMin > epsilon &&
                particle.Position.X < body.MinPosition.X &&
                body.MinPosition.Z - r2 < particle.Position.Z && particle.Position.Z < body.MaxPosition.Z + r2 &&
                body.MinPosition.Y - r2 < particle.Position.Y && particle.Position.Y < body.MaxPosition.Y + r2)
            {
                norm = new Vector3(-1, 0, 0);
                adjustment = stiff * diffXMin - damp * DotProduct(norm, particle.VelocityEvaluative);
                Vector3 acc = norm * adjustment;
                particle.Acceleration += acc;
            }

            if (diffXMax > epsilon &&
                body.MaxPosition.X < particle.Position.X &&
                body.MinPosition.Z - r2 < particle.Position.Z && particle.Position.Z < body.MaxPosition.Z + r2 &&
                body.MinPosition.Y - r2 < particle.Position.Y && particle.Position.Y < body.MaxPosition.Y + r2)
            {
                norm = new Vector3(1, 0, 0);
                adjustment = stiff * diffXMax - damp * DotProduct(norm, particle.VelocityEvaluative);
                Vector3 acc = norm * adjustment;
                particle.Acceleration += acc;
            }

            if (diffZMin > epsilon &&
                particle.Position.Z < body.MinPosition.Z &&
                body.MinPosition.X - r2 < particle.Position.X && particle.Position.X < body.MaxPosition.X + r2 &&
                body.MinPosition.Y - r2 < particle.Position.Y && particle.Position.Y < body.MaxPosition.Y + r2)
            {
                norm = new Vector3(0, 0, -1);
                adjustment = stiff * diffZMin - damp * DotProduct(norm, particle.VelocityEvaluative);
                Vector3 acc = norm * adjustment;
                particle.Acceleration += acc;
            }

            if (diffZMax > epsilon &&
                body.MaxPosition.Z < particle.Position.Z &&
                body.MinPosition.X - r2 < particle.Position.X && particle.Position.X < body.MaxPosition.X + r2 &&
                body.MinPosition.Y - r2 < particle.Position.Y && particle.Position.Y < body.MaxPosition.Y + r2)
            {
                norm = new Vector3(0, 0, 1);
                adjustment = stiff * diffZMax - damp * DotProduct(norm, particle.VelocityEvaluative);
                Vector3 acc = norm * adjustment;
                particle.Acceleration += acc;
            }


            //diff = 2 * radius - (particle.Position.Y - body.MaxPosition.Y) * ss;
            //if (diff > epsilon)
            //{
            //    norm = new Vector3(0, 1, 0);
            //    adjustment = stiff * diff - damp * DotProduct(norm, particle.VelocityEvaluative);
            //    acceleration += norm * adjustment;
            //}

            //diff = 2 * radius - (body.MinPosition.Y - particle.Position.Y) * ss;
            //if (diff > epsilon)
            //{
            //    norm = new Vector3(0, -1, 0);
            //    adjustment = stiff * diff - damp * DotProduct(norm, particle.VelocityEvaluative);
            //    acceleration += norm * adjustment;
            //}

        }

        float DotProduct(Vector3 left, Vector3 right)
        {
            return left.X * right.X +
                   left.Y * right.Y +
                   left.Z * right.Z;
        }

        private void SPH_ComputePressureGrid(int particleIndex)
        {

            //Parallel.For(0, particles.Count, i =>
            //for (int i = 0; i < particles.Count; i++)
            {
                float dx, dy, dz, sum, dsq, c;
                float d, d2, mR, mR2;
                float radius = Parameters.SmoothingRadius / Parameters.SimScale;
                d = Parameters.SimScale;
                d2 = d * d;
                mR = Parameters.SmoothingRadius;
                mR2 = mR * mR;
                var particle = particles[particleIndex];
                sum = 0;
                if (particle.IsFixed)
                    return;

                neighborCount[particleIndex] = 0;
                ushort count = 0;
                var gridCells = Grid_FindCells(particle.Position, radius);

                foreach (int cellValue in gridCells)
                {
                    if (cellValue <= -1 || cellValue >= grid.Count)
                    {
                        continue;
                    }
                    var otherParticleIndex = grid[cellValue];
                    while (otherParticleIndex != -1)
                    {
                        var otherParticle = particles[otherParticleIndex];
                        if (particle == otherParticle)
                        {
                            otherParticleIndex = otherParticle.Next;
                            continue;
                        }
                        //distance in cm
                        dx = (particle.Position.X - otherParticle.Position.X) * d;
                        dy = (particle.Position.Y - otherParticle.Position.Y) * d;
                        dz = (particle.Position.Z - otherParticle.Position.Z) * d;
                        dsq = (dx * dx + dy * dy + dz * dz);
                        if (mR2 > dsq)
                        {
                            //h^2 - r^2
                            c = mR2 - dsq;
                            //(h^2 - r^2)^3
                            sum += c * c * c;

                            //Save neighbor particle for later calculations
                            if (count < maxNeighbors)
                            {
                                neighbor[particleIndex][count] = (ushort)otherParticleIndex;
                                neighborDist[particleIndex][count] = (float)Math.Sqrt(dsq);
                                count++;
                                neighborCount[particleIndex] = count;
                            }
                        }
                        otherParticleIndex = otherParticle.Next;
                    }

                } //end looping over cell
                particle.Density = sum * Parameters.ParticleMass * poly6Kern;
                particle.Pressure = (particle.Density - Parameters.RestDensity) * Parameters.InternalStiffness;

                particle.Density = Parameters.RestDensity;
                //particle.Density = 1.0f / particle.Density;
            } // end looping over particles
            //);
        }

        private void SPH_ComputeForceGridNC(int particleIndex)
        {
            //Parallel.For(0, particles.Count, i =>
            //for (int i = 0; i < particles.Count; i++)
            {
                Vector3 force;
                float pterm, vterm, dterm;
                float c, d;
                float dx, dy, dz;
                float mR, mR2, visc;

                d = Parameters.SimScale;
                mR = Parameters.SmoothingRadius;
                mR2 = mR * mR;
                visc = Parameters.Viscosity;
                var particle = particles[particleIndex];
                if (particle.IsFixed)
                    return;

                force = new Vector3();

                for (int j = 0; j < neighborCount[particleIndex]; j++)
                {
                    var otherParticle = particles[neighbor[particleIndex][j]];
                    float otherPressure = otherParticle.Pressure;
                    float otherDensity = otherParticle.Density;
                    Vector3 otherVelocityEvaluative = otherParticle.VelocityEvaluative;
                    if (otherParticle.IsFixed)
                    {
                        //otherDensity = particle.Density;
                        //otherPressure = particle.Pressure;
                        //otherVelocityEvaluative = particle.VelocityEvaluative;
                    }
                    //distance in cm
                    dx = (particle.Position.X - otherParticle.Position.X) * d;
                    dy = (particle.Position.Y - otherParticle.Position.Y) * d;
                    dz = (particle.Position.Z - otherParticle.Position.Z) * d;

                    float dist = neighborDist[particleIndex][j];
                    c = (mR - dist);

                    pterm = 0.5f * c * spikyKern * (particle.Pressure + otherPressure) / dist;
                    if (particle.Density == 0 || otherDensity == 0)
                        dterm = 0;
                    else
                        dterm = c / ( particle.Density * otherParticle.Density);
                    vterm = laplacianKernel * visc;
                    force.X += (pterm * dx + vterm * (otherVelocityEvaluative.X - particle.VelocityEvaluative.X)) * dterm;
                    force.Y += (pterm * dy + vterm * (otherVelocityEvaluative.Y - particle.VelocityEvaluative.Y)) * dterm;
                    force.Z += (pterm * dz + vterm * (otherVelocityEvaluative.Z - particle.VelocityEvaluative.Z)) * dterm;

                } // end looping neighbor particles
                particle.SPHForce = force;
            } //end looping particles
            //);
        }

        private void ComputeKernels()
        {
            Parameters.ParticleDistance = (float)Math.Pow(Parameters.ParticleMass / Parameters.RestDensity, 1 / 3.0);
            smoothing2 = Parameters.SmoothingRadius * Parameters.SmoothingRadius;

            // Wpoly6 kernel (denominator part) - 2003 Mueller, p.4
            poly6Kern = (float)(315.0 / (64.0 * Math.PI * Math.Pow(Parameters.SmoothingRadius, 9)));

            // Laplacian of viscocity (denominator): PI h^6
            spikyKern = (float)(45.0 / (Math.PI * Math.Pow(Parameters.SmoothingRadius, 6)));

            laplacianKernel = (float)(45.0 / (Math.PI * Math.Pow(Parameters.SmoothingRadius, 6)));
        }

        public override void Reset()
        {
            particles.Clear();
            DeltaTime = 0.003f;
            DeltaTime = 0.0015f;
            if (this.Parameters == null)
                this.Parameters = GetDefaultParameters();
            //base.Reset();

            this.Parameters.PlaneGravityEnabled = true;
            this.Parameters.PointGravityEnabled = false;

            this.Parameters.PointGravityPosition = new Vector3(0, 0, 50);
            this.Parameters.PlaneGravityDirection = new Vector3(0, 0, -9.8f);

        }

        public override int AddPoint()
        {
            FluidParticle particle = new FluidParticle();
            int index = particles.Count;

            particle.SPHForce = new Vector3();
            particle.Velocity = new Vector3();
            particle.VelocityEvaluative = new Vector3();
            particle.Next = 0;
            particle.Pressure = 0;
            particle.Density = 0;

            particles.Add(particle);

            return index;
        }

        public override int AddPointReuse()
        {
            FluidParticle particle;
            int index = particles.Count;
            if (particles.Count < Parameters.NumParticles)
            {

                particle = new FluidParticle();
                particles.Add(particle);
            }
            else
            {
                index = random.Next(0, index - 1);
                particle = particles[index];
            }

            particle.SPHForce = new Vector3();
            particle.Velocity = new Vector3();
            particle.VelocityEvaluative = new Vector3();
            particle.Next = 0;
            particle.Pressure = 0;
            particle.Density = 0;

            return index;
        }

        //------------------------------------------------------ SPH Setup 
        //
        //  Range = +/- 10.0 * 0.006 (r) =	   0.12			m (= 120 mm = 4.7 inch)
        //  Container Volume (Vc) =			   0.001728		m^3
        //  Rest Density (D) =				1000.0			kg / m^3
        //  Particle Mass (Pm) =			   0.00020543	kg						(mass = vol * density)
        //  Number of Particles (N) =		4000.0
        //  Water Mass (M) =				   0.821		kg (= 821 grams)
        //  Water Volume (V) =				   0.000821     m^3 (= 3.4 cups, .21 gals)
        //  Smoothing Radius (R) =             0.02			m (= 20 mm = ~3/4 inch)
        //  Particle Radius (Pr) =			   0.00366		m (= 4 mm  = ~1/8 inch)
        //  Particle Volume (Pv) =			   2.054e-7		m^3	(= .268 milliliters)
        //  Rest Distance (Pd) =			   0.0059		m
        //
        //  Given: D, Pm, N
        //    Pv = Pm / D			0.00020543 kg / 1000 kg/m^3 = 2.054e-7 m^3	
        //    Pv = 4/3*pi*Pr^3    cuberoot( 2.054e-7 m^3 * 3/(4pi) ) = 0.00366 m
        //     M = Pm * N			0.00020543 kg * 4000.0 = 0.821 kg		
        //     V =  M / D              0.821 kg / 1000 kg/m^3 = 0.000821 m^3
        //     V = Pv * N			 2.054e-7 m^3 * 4000 = 0.000821 m^3
        //    Pd = cuberoot(Pm/D)    cuberoot(0.00020543/1000) = 0.0059 m 
        //
        // Ideal grid cell size (gs) = 2 * smoothing radius = 0.02*2 = 0.04
        // Ideal domain size = k*gs/d = k*0.02*2/0.005 = k*8 = {8, 16, 24, 32, 40, 48, ..}
        //    (k = number of cells, gs = cell size, d = simulation scale)
        public FluidSystemParams GetDefaultParameters()
        {
            FluidSystemParams parameters = new FluidSystemParams()
            {
                NumParticles = 1024,
                SimScale = 0.010f,
                Viscosity = 0.2f,
                RestDensity = 600.0f,
                ParticleMass = 0.00020543f,
                ParticleRadius = 0.004f,
                ParticleDistance = 0.0059f,
                SmoothingRadius = 0.01f,
                InternalStiffness = 0.5f,
                ExternalStiffness = 20000.0f,
                ExternalDamping = 256.0f,
                Limit = 200.0f,
                PlaneGravityDirection = new Vector3(0, 0, -9.8f),
                PlaneGravityEnabled = true,
                PointGravityEnabled = false
            };

            parameters.InternalStiffness = 2f;
            parameters.RestDensity = 200f;
            //parameters.SmoothingRadius = 0.005f;
            //parameters.ParticleMass *= 1.0f;
            //parameters.Viscosity = 0.4f;
            parameters.ParticleDistance *= 1.25f;
            parameters.ParticleRadius *= 1.250f;
            return parameters;
        }

        #endregion
    }
}
